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Abstract 

We propose rotation inferred from the polar decomposition of the 
flow gradient as a diagnostic for elliptic (or vortex-type) invariant re¬ 
gions in non-autonomous dynamical systems. We consider here two- 
and three-dimensional systems, in which polar rotation can be char¬ 
acterized by a single angle. Eor this polar rotation angle (PRA), we 
derive explicit formulas using the singular values and vectors of the 
flow gradient. We find that closed level sets of the PRA reveal ellip¬ 
tic islands in great detail, and singular level sets of the PRA uncover 
centers of such islands. Both features turn out to be objective (frame- 
invariant) for two-dimensional systems. We illustrate the diagnostic 
power of PRA for elliptic structures on several examples. 


1 Introduction 

Complex dynamical systems exhibit a mixture of chaotic and coherent be¬ 
havior in their phase space. The latter manifests itself in coherent islands of 
regular behavior surrounded by a chaotic background flow. The best known 
classic examples of such islands are formed by Kolmogorov-Arnold-Moser 
(KAM) tori, composed of quasi-periodic trajectories in Hamiltonian systems 
[see, e.g., 1,2]. Outside elliptic regions filled by such tori, chaotic trajectories 
dominate the dynamics. 
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Even more intriguing is the existence of similar elliptic islands in tur¬ 
bulent fluid flow, as broadly confirmed by experiments and numerical sim¬ 
ulations [see, e.g., 3, 4]. Just as KAM islands, coherent vortices capture 
trajectories and keep them out of chaotic mixing zones. Unlike KAM tori, 
however, coherent vortices are composed of trajectories that are generally 
not recurrent in any frame. During their finite time of existence, these co¬ 
herent vortices traverse without filamentation but also without displaying 
any particular periodic or quasiperiodic pattern. Still, we generally refer 
to such regions here as elliptic, as they mimic the dynamic role of elliptic 
islands occupied by classic KAM tori. 

Eulerian approaches to describing elliptic islands seek domains where 
rotation dominates the instantaneous velocity field. At the simplest level, 
this involves locating regions of closed streamlines, high enough vorticity or 
low enough pressure (cf. Jeong and Hussain [5] and Haller [6] for reviews). 
Such domains reveal instantaneous velocity field features at a low cost, but 
are unable to frame long-term material coherence exhibited by trajectories. 
In addition, the results from these instantaneous approaches depend on the 
choice of scalar thresholds and on the frame of reference. 

More sophisticated Eulerian principles for elliptic regions seek sets of 
points where rotation dominates strain (see, e.g., Jeong and Hussain [5], 
Okubo [7], Weiss [8], Hunt et al. [9], Hua and Klein [10], Tabor and Klap- 
per [11], and also Jeong and Hussain [5] and Haller [6] for reviews). These 
principles infer both rotation and strain from the instantaneous velocity 
gradient, thereby rendering the results Galilean invariant. The elliptic re¬ 
gions they provide, however, still change under rotations of the frame. Since 
truly unsteady flows have no distinguished frame of reference [12], frame- 
dependence in the detection of vortical structures is an impediment. Indeed, 
the available measurement velocity data of geophysical flows is often given 
in a rotating frame to begin with, and no optimal frame is known a priori 
for structure detection. More importantly, no mathematical relationship is 
known (or likely to exist) between instantaneous rotation-strain principles 
and material coherence over extended time intervals. 

In contrast, Lagrangian approaches to elliptic islands seek to identify re¬ 
gions where trajectories stay close for longer periods. These approaches can 
roughly be divided into three categories: geometric, set-based and diagnos¬ 
tic methods. The geometric methods identify elliptic domain boundaries as 
spacial closed material lines showing no filamentation [13-15] or curvature 
change [16]. Set-based methods partition the phase space into almost in¬ 
variant subsets (see Budisic et al. [17], Eroyland [18] and references therein). 
While the boundaries of such sets may undergo filamentation, the over- 
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all subsets remain largely coherent. Finally, diagnostic approaches propose 
Lagrangian scalar fields whose features are expected to distinguish mixing 
regions from coherent ones [19-24]. These Lagrangian methods do not re¬ 
turn identical results and are not backed by specific mathematical results on 
the features they highlight. In fact, the material invariance of the extracted 
vortical boundaries is only guaranteed in the case of the geodesic approach 
of Haller and Beron-Vera [13] and Haller [15]. 

The Lagrangian methods listed above focus on stretching or lack thereof. 
In contrast, very few Lagrangian diagnostics target rotation, even though 
sustained and coherent rotation is perhaps the most striking feature of tra¬ 
jectories forming elliptic islands. One of the few exceptions targeting mate¬ 
rial rotation is the finite-time rotation number (FTRN), developed to detect 
hyperbolic (i.e., repelling or attracting as opposed to vortical) structures 
through its ridges [25]. The FTRN assumes that the dynamical system is 
defined via an iterated map with an annular phase space. For dynamical sys¬ 
tems with general time dependence and non-annular phase space, however, 
this approach is not applicable. This also means that the approach is frame- 
dependent, given that translations and rotations will generally destroy the 
time-periodicity of a dynamical system. 

Another Lagrangian diagnostic involving a consideration of rotation is 
the mesocronic analysis of Mezic et al. [23]. This approach offers a formal 
extension of the Okubo-Weiss principle from the velocity gradient to the 
flow gradient, classifying an initial condition as elliptic if the flow gradi¬ 
ent has complex eigenvalues at that point. The mesoelliptic diagnostic is 
efficient to compute and has been shown to mark vortical regions in sev¬ 
eral cases. The direct extension from the Okubo-Weiss principle, however, 
also renders the mesoelliptic diagnostic frame-dependent. In addition, the 
complex eigenvalues of a finite-time flow map have no known mathemati¬ 
cal relationship with elliptic islands in flows with general time dependence. 
Accordingly, some annular subsets of classic elliptic domains fail the test of 
meso-ellipticity even in steady flows (cf. [23], Fig. 1). 

Here we propose a mathematically precise assessment of material rota¬ 
tion, the polar rotation angle (PRA), as a new diagnostic for elliptic islands 
in two- and three-dimensional flows. The PRA is the angle of the rigid- 
body rotation component obtained from the classic polar decomposition of 
the flow gradient into a rotational and a stretching factor. We show how the 
PRA can readily be computed from invariants of the flow gradient and the 
Cauchy-Green strain tensor. Level sets of the PRA turn out to be objective 
(frame-invariant) in planar flows. We find that these level sets reveal the 
internal structure of elliptic islands in great detail at a relatively low com- 
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putational cost. We also find that local extrema of the PRA mark elliptic 
island centers suitable for automated vortex tracking in Lagrangian fluid 
dynamics. 


2 Preliminaries 


2.1 Set-up 

Consider the dynamical system 


X = u(x, t), X G P C t G / C M, 


( 1 ) 


with the corresponding flow map 



Xo H-). x(t;to,xo), 


( 2 ) 


the diffeomorphism that takes the initial condition xq to its time-t position 
x(t;to,xo) under system (1). Here, V denotes the phase space and / is a 
finite time interval of interest. 

The deformation gradient governs the infinitesimal deformations 

of the phase space V. In particular, an initial perturbation ^ at point xq 
and time to is mapped, under the system (1), to VF^q(xo)^ at time t. We 
also define the Cauchy-Green strain tensor^ 




(3) 


where the symbol T denotes matrix transposition. The tensor C^q(xo) is 
symmetric and positive definite. Therefore, it has an orthonormal set of 
eigenvectors {^i(xo),^ 2 (xo),^ 3 (xo)}. The corresponding eigenvalues 0 < 
Ai(xo) < A 2 (xo) < A 3 (xo) therefore satisfy 


C^xo)^j(xo) = Aj(xo)^i(xo), i E {1,2,3}, 

(^i(xo),^fc(xo)) = 0, j, fee (1,2,3}, jT^k, 


(4) 

(5) 


with (•, •) denoting the Euclidean inner product. For notational simplicity, 
we omit the dependence of the eigenvalues and eigenvectors on to and t. Also, 
we consider two-dimensional flows as a special case satisfying dx^Ui{'K^ t) = 0, 
i = 1,2,3. 
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2.2 Polar decomposition 

Any square matrix admits a factorization into the product of a unitary ma¬ 
trix with a symmetric positive-semidefinite matrix [26]. When the square 
matrix is nonsingular, such as VF^^, then the symmetric factor in the de¬ 
composition is positive definite. 

Specifically, the deformation gradient admits a unique decomposi¬ 

tion of the form 

= ( 6 ) 

where the 3x3 matrices and have the following properties [26-28]: 

1. The rotation tensor is proper orthogonal, i.e., 

(RUX=R.‘.«)^=I. detR‘. = l. 

2. The right stretch tensor is symmetric and positive-definite, satis- 
fying 

[uy"=c‘„, (7) 

3. The eigenvalues of are with corresponding eigenvectors 

U^xo)Cfc(xo) = A/Afc(xo)^fc(xo), k = 1,2,3, (8) 


4. The time derivative of the rotation tensor satisfies 


R!.= (w(x(t),i)--R‘„ u*„(u‘„)-'-(u‘„) u‘„ 

( 9 ) 


where W = ^ 


Vu-(Vu 




is the vorticity (or spin) tensor and x(t) 

is a shorthand notation for the trajectory x(t;to,xo). A derivation of 
(9) can be found, e.g., in [29, Section 23]. 


The geometric interpretation of the polar decomposition is the following 
[27, 30]. At any point xq of the phase space, the orthogonal basis {^k}i<k<3 
is mapped into {^Fl^{-xo)^k}i<k<3 under the linearized flow map VF^^. Any 
stretching and compression in the deformation is encoded into the stretch 
tensor while the overall rigid-body rotation of material elements is 

encoded into the rotation tensor R^^. Figure 1 illustrates the action of these 
tensors on area elements in two dimensions. 
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Figure 1: The action of the deformation gradient is uniquely decom¬ 

posable into positive definite stretch by followed by rotation by 
This results in the polar decomposition 


When mapped forward under the deformation gradient VF^^, a general 
unit vector a experiences two stages of rotation. First, a is rotated (and 
simultaneously stretched) by the stretch tensor into the vector U^^a. 
This first stage of rotation is entirely due to shear, with the magnitude and 
axis of rotation depending on a. The second stage of rotation experienced 
by a is due to the rotation tensor R^^, which rotates a into its final position 
R^^U^^a at time t. This second rotation acts in the same way on all U^^a 
vectors by the proper orthogonal nature of R^^. 

Formed by the eigenvectors of the principle rectangle illustrated 
in Fig. 2 has a special feature: it is the unique rectangle on which the 
first stage of rotation under is inactive. This is because the edges of the 
principal rectangle align with the eigenvectors of (cf. eq. (8) and Fig. 1), 
and hence remain unrotated by The total rotation experienced by the 
edges of the principal rectangle is, therefore, just the rigid-body rotation 
exerted by the rotation tensor R^^. 

Formula (9) shows the difference between instantaneous Eulerian rota¬ 
tion measured by the vorticity tensor W and finite-time material rotation 
measured by R^^. In particular, at an initial time to, we have 

^toL=to 

but differs from the vorticity tensor W(x(t),t) for times t ^ to- 
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VF‘(x„)?,(x„) 



Figure 2: Finite-time deformation of an area element of a two-dimensional 
phase space under the flow map The orthonormal basis {^ 1 ,^ 2 } is 

mapped to the orthogonal basis Other initially orthogo¬ 

nal material elements, such as the diagonals shown in blue, are mapped to 
non-orthogonal material elements. 

3 Polar rotation angle (PRA) 

The classic procedure for computing the polar decomposition in continuum 
mechanics starts with determining U^^as the principal square root of the 
Cauchy-Green strain tensor (cf. formula (7)). This is the simplest to do by 
diagonalizing taking the positive square root of its diagonal elements, 
and transforming back the resulting matrix from the strain eigenbasis to the 
original basis. Next, one obtains the rotation tensor directly from (6) as 
[U^o] More efficient numerical procedures are also available 
(see [31] and the references cited therein) 

These computational approaches, however, offer little insight into the ge¬ 
ometry of the rotation generated by . Taking a more geometric approach, 
one may recall that any three-dimensional rotation R^^ has a Rodrigues rep- 


7 




resentation [32] of the form 

= I + + (1 - cosel) [Py ^ (10) 

where I is the 3x3 identity matrix and is a 3 x 3 skew-symmetric matrix 
such that 

The unit vector is the eigenvector of corresponding to its unit eigen¬ 
value, i.e., 

P'Uxo)ryxo) = ryxo). (11) 

For planar flows, the eigenvector is the unit normal to the plane of 
motion, and hence is independent of xq. In three-dimensions, depends 
on the location xq in a way discussed in the next section (cf. Proposition 1) 
Once an orientation for the unit vector is selected, the angle 

0^q(xo) G [0,27r) is uniquely determined. This angle measures the amount 
of local solid-body rotation experienced by material elements along the tra¬ 
jectory x(t;to,xo). 

Definition 1. We refer to the scalar function 

<9*0 (xo) G [0,27r) 

determined by (10) as the polar rotation angle (PRA) at the initial condition 
Xq with respect to the time interval 


4 Computing the PRA 


Taking the trace of both sides in (10), then taking the skew-symmetric part 
of both sides of (10) yields the formulas 


cos 


= O - 1) ’ 




K 


Mj 


[Pt] 







(12a) 

(12b) 


To evaluate the expression for cos^^^ in (12a), Guan-Suo [33] expressed 
tr as a somewhat cryptic function of the scalar invariants of the matrices 
\ ^to' we derive a simply computable 

and intuitive alternative that only involves quantities arising in typical La- 
grangian coherent structure calculations [15]: the deformation gradient, and 
the eigenvalues and eigenvectors of the the Cauchy-Green strain tensor. 





Proposition 1. 

(1) In three-dimensional flows, the PRA satisfies the relations 


= - 






sin 0*0 = 




(13a) 

i/jG {1,2,3}, (13b) 


where e = ( 61 , 62 , 63 )^ is the normalized eigenvector corresponding to 
the unit eigenvalue of the matrix 





j,k G {1,2,3}, 


and Cijk is the Levi-Civita symbol. Furthermore, we have e = where 
is the axis of rotation defined by ( 11 ). 

(2) In two-dimensional flows, we have 


^ z = 1 or 2, (14a) 

V 

sin 0*0 ^ (- 1 )^ (i,j) = (l, 2 ) or ( 2 , 1 ), 

W 

where Ai < A 2 are the eigenvalues of the two-dimensional Cauchy- 
Green strain tensor with corresponding eigenvectors and ^ 2 - 



Proof. See Appendix A. □ 

Using both expressions in the formulas (13) (or formulas (14), in the 
two-dimensional case), the four-quadrant polar rotation angle G [0, 27r) 
can be reconstructed as 

^to = [1 - sign (sin0**J] tt + sign (sin 04 *J cos“^ (cos0tj , (15) 


where 


sign((a) 


1 

-1 


if (a > 0 
if (a < 0 


For completeness, in Appendix B, we also derive a formula for the total 
rotation of an arbitrary material element, not just for the strain eigenvectors. 
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Evaluating this general formula is computationally more costly, as it involves 
advecting initial directions by the flow map through all intermediate times 
within the interval In addition, due to the non-rigid-body nature of 

deformation along a trajectory, the total material rotation will be different 
for different material elements. When evaluated on initial directions aligned 
with and ^2, however, this total Lagrangian rotation agrees with the PRA 
modulo multiples of 27t. 

5 Polar LCS 

A recent approach to the systematic detection of elliptic Lagrangian coherent 
structures (LCS) targets closed material lines that exhibit no filamentation 
over the finite time interval [to,t] (Haller and Beron-Vera [13], Haller [15]). 
These elliptic LCSs turn out to be uniformly stretching closed material lines, 
i.e., all their subsets exhibit the same relative stretching. Outermost mem¬ 
bers of nested elliptic LCS families then serve as the ideal boundaries of 
perfectly coherent elliptic islands. 

Here we propose a dual approach to elliptic LCSs by requiring uniformity 
in the polar rotation of material elements forming the LCS, as opposed to 
uniformity in their stretching. 

Definition 2. A polar Lagrangian coherent structure (polar LCS) over the 
time interval [to,t] is a closed (i.e., tubular in 3D and circular in 2D) and 
connected codimension-one material surface whose time to position is a level 
set of (9^X0). 

As any material surface, a polar LCS is invariant under the flow. It 
is formed by trajectories starting from a closed and connected level set of 
0 ^q(xo) at time to. The following simple observation shows that polar LCSs 
can be detected as connected and closed level sets of trigonometric functions 
of 0 ^q(xo), and hence are directly computable from the formulas (13)-(14). 

Proposition 2. Connected components of the level sets of cos 6l^ and 
sin coincide with connected components of the level sets of 0 |q(xo). 

Proof. Assume the contrary, i.e., the existence of two points xq and xq that 
are in the same connected component of a level set C of cos9l^{'Xo) but on 
different connected level sets of ^^^(xo). Then on any continuous path con¬ 
necting Xq and Xq, the polar rotation angle 9l^ should change continuously 
from 0 |q(xo) to ^^^(xo) ^ 0 |q(xo), and hence cos^^^ cannot be constant along 
this path. Since xq and xq are in the connected set £, there is therefore a 
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continuous path connecting xq and xq within C along which cos0^q(xo) can¬ 
not be constant. But this contradicts the assumption that £ is a level set 
of cos0^q(xo). The argument for sin^^^ is identical. □ 

A practical consequence of Proposition 2 is that connected level sets of 
0 ^q(xo) can be constructed as those of cos^l^ and sin^^^, without verifying 
the orientability of on U. This renders the computation of the tensor 

and the rotation axis unnecessary, as one can compute the two- 
quadrant angle 6l^ G [0,7r] from equation (13a) as 


el = cos ^ 






(16) 


Proposition 2 ensures that the level sets of 9\^ computed from (16) coincide 
with those of the four-quadrant PRA angle computed from (15). 

All quantities derived from the deformation gradient are invariant 

with respect to time-dependent translations of the coordinate frame. There¬ 
fore, polar LCSs are Galilean invariant objects. For two-dimensional flows, 
polar LCSs also turn out to be invariant under time-dependent rotations of 
the frame. In the language of continuum mechanics [29], polar LCSs in two 
dimensions are objective. 

Proposition 3. In two-dimensional flows, a polar LCS over the time interval 
[to,t] is objective, i.e., invariant under coordinate changes of the form 


X = Q(i)y+ b(f), (17) 

where Q(t) € 50(2) and b(f) € are smooth functions of time t. 

Proof. See Appendix C. □ 

An elliptic island marked by the PRA has a natural center point: the 
PRA extremum point surrounded by closed PRA contours. This leads to 
the following definition of a Lagrangian vortex center: 

Definition 3. A Lagrangian vortex center over the time interval [to^t] is a 
set of trajectories evolving from a connected, codimension-two level set of 

A Lagrangian vortex center identified from the PRA is, therefore, com¬ 
posed of a single trajectory in two dimensions and of a one-parameter fam¬ 
ily of trajectories (i.e., a material line) in three dimensions. Despite recent 
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progress in the accurate detection of coherent Lagrangian vortex boundaries 
[15], approaches to Lagrangian vortex center definition and detection have 
notably been missing. As we illustrate in Section 6.2 below, vortex centers 
defined as PRA extrema indeed show distinguished behavior: they capture 
the translational motion of an elliptic island without being affected by the 
rotational motion of trajectories inside the island. As connected level sets 
of the PRA, the Lagrangian vortex centers defined in Definition 3 are also 
objective in two-dimensional flows (cf. Proposition 3). 

6 Examples 

In this section, we compute the PRA on several examples to illustrate how 
its closed level curves (i.e., initial positions of polar LCSs) highlight the 
internal structure of elliptic islands in detail. 

6.1 Standard map 

We first consider the standard map 

In+l ~ ^ sin 072, 

0n+l ~ 4^n “1“ (1^) 

which is a Poincare map 'P of a rotor excited by a periodic impulsive 
force [34]. In the absence of the impulse, i.e., for e = 0, the angular mo¬ 
mentum In is constant and the angular position (^n increases linearly as an 
integer multiple of the angular momentum. 

For 6^0, however, the system can exhibit complicated dynamics. De¬ 
pending on the initial condition (/o, 0 o )5 the trajectories may be periodic, 
quasi-periodic or chaotic. The quasi-periodic trajectories lie on KAM tori, 
the classic examples of vortical structures that we wish to visualize through 
the PRA. 

The left plot in Fig. 3a shows 1000 iterations of the standard map for 400 
uniformly distributed initial conditions and e = 1. This reveals invariant 
KAM tori, resonance islands and chaotic regions. The right panel of the 
same figure shows the PRA, computed from formula (15) with i = 2, with 
the flow map being equal to 200 iterations of the map (18), i.e. (xo) = 
P^°°(xo), where xq = (^oVo) and ((/>n+i,/n+i) = 'P{4’n,In)- To ensure the 
accuracy of the finite differences for the computation of the deformation 
gradient VF^^, we use a dense grid of initial conditions consisting of 1000 x 
1000 uniformly distributed points over the phase space = [0, 27r] x [0, 27r]. 
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Figure 3: Standard map. (a) Left: 1000 iterations of the standard map for 
400 uniformly distributed initial conditions over the torus [0, 27r] x [0,27r]. 
Right: The PRA 9l^ for 200 iterations of the standard map, clearly marking 
polar LCSs (closed contours) and Lagrangian vortex centers (local extrema) 
(b) The close-up view of the region marked by a rectangle in (a). 


Figure 3b shows a close-up of a region of the phase space containing 
chaotic trajectories, KAM tori and a period-5 resonance island. For this 
close-up view, the Poincare map is recomputed from 1000 iterations of 2500 
initial conditions. The corresponding PRA plot on the right is computed 
only from 500 iterations, i.e., from 

We conclude that the KAM tori and resonance islands are sharply en¬ 
hanced by the PRA relative to a simple iteration of the map, even though 
the number of iterations used in constructing the PRA plot is only half the 
number used for the Poincare map. The chaotic region is marked by small- 
scale rapid variations in the PRA, in line with the sensitive dependence of 
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X 

Figure 4: Left: Vorticity uo at the initial time t — 50. Right: The PRA 9\^ 
computed from formula (15) over the time interval [50,100]. 



the rotation angle on initial conditions in these regions. 

Figure 3 also shows that the center-type fixed points in the elliptic islands 
are clearly marked with local extrema of the PRA, supporting the idea of 
defining Lagrangian elliptic centers as stated in Definition 3. 


6.2 Two-dimensional turbulence 

Consider the Navier-Stokes equation 

dtw + u • Vu = -Vp + z/Au + f, V • u = 0, (19) 


where z/ is the kinematic viscosity and f denotes the forcing. For an ideal 
two-dimensional fluid flow (z/ = 0 and f = 0 ), the vorticity cj, given by 
V X u = (0,0,cj), is preserved along fluid trajectories, i.e.. 


Dcj 


( 20 ) 


where ^ + u • V is the material derivative. Therefore, closed level 

curves of vorticity are material curves, acting as barriers to the transport of 
fluid particles. In the presence of molecular diffusion and external forcing, 
however, vorticity is not a material invariant and hence its closed contours 
no longer signal elliptic islands for fluid trajectories. 

To illustrate the use of PRA in detecting elliptic islands in a turbulent 
flow, we solve the Navier-Stokes equation (19) with v = 10“^ on the domain 
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V = [0, 27 r] X [0, 27 r] with periodic boundary conditions. We use a pseudo- 
spectral method with 512^ modes to evaluate the spatial partial derivatives 
and the nonlinear term. The external forcing is random in phase and only 
active over the wave-numbers 3.5 < k < 4.5. The forcing amplitude is time- 
dependent and chosen to balance the instantaneous enstrophy dissipation 
—z/ |V6c;(x, t)|^dx. The time integration is carried out by a variable step- 
size, fourth-order Runge-Kutta method [35]. We solve the equation up to 
time t = 100. We observe that the turbulent flow is fully developed after 50 
time units. Therefore, we choose times to = 50 and t = 100 as the initial 
and final times for the computation of the PRA 9l ^. 

Such two-dimensional turbulent flows tend to generate long-lasting co¬ 
herent vortices [36], which are also prevalent in geophysical flows [37]. Highly 
coherent Lagrangian signatures of such vortices have been recently identified 
as regions bounded by uniformly stretching material lines [13, 15]. 

Here, we take an alternative approach and identify coherent Lagrangian 
vortices as regions filled with polar LCSs. In other words, we seek the elliptic 
islands of turbulence as regions of closed material lines that pointwise have 
the same rigid-body rotation component in their deformation over the time 
interval of interest. 

Figure 4 (right panel) shows the PRA computed from formula (15) for 
512 X 512 uniformly distributed initial conditions. The polar LCSs are clearly 
visible as concentric closed contours of ^^^(xo). Figure 5 shows a closeup 
view of a coherent Lagrangian vortex identified from the PRA plot. Note 
how the PRA shows a sharp distinction between the vortical region and the 
surrounding chaotic background. As in the case of the standard map (see 
Fig. 3), the chaotic region is marked by small-scale, sharp variations of the 
Lagrangian rotation due to sensitive dependence of material rotation angle 
on initial conditions. 

While the velocity field is well-resolved, resolving small-scale Lagrangian 
structures requires significantly higher resolution [38-40]. At the present 
resolution, the Lagrangian structures in the chaotic region are not well- 
resolved. Nonetheless, the boundary of the vortex can be approximated 
by the contour across which the PRA transitions from concentric large- 
scale contours to small-scale sharp variations (see the red-colored contour in 
Fig. 5). 

We now illustrate that the large-scale polar LCSs, defined by closed con¬ 
tours of the PRA (cf. Definition 2) indeed remain coherent under advection. 
We advect two such contours under the flow, with their advected positions 
shown in the right panel on Fig. 5 at time t = 100. 

As a measure of coherence we define relative stretching of material lines 
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Figure 5: Left: Contours of the PRA signaling coherent and chaotic regions. 
Right: Advected image of select contours to the final time t — 100. The local 
extremum of the PRA (marked by a cross) defining the Lagragian vortex 
center by Definition 3. 

as [l{t) — ^(to)] /^(^o)? where I denotes the length of the material line as a 
function of time. The relative stretching of the blue and red contours are 
2.65% and —1.38%, respectively. These relative stretching values remain in 
the order of stretching exhibited by perfectly coherent elliptic LCSs obtained 
from the geodesic LCS theory [13]. 




Figure 6: Left panel: Trajectories of the Lagrangian vortex center (red) 
and nearby passive tracers (blue and black). Middle and Right panels: The 
coordinates of the vortex center and nearby tracers as a function of time. 

The cross in Fig. 5 marks a local extremum of the PRA, which is a 
Lagrangian vortex center by our Definition 3. This local extremum indeed 
turns out to behave as the vortex center over the time interval of interest, 
i.e., t G [50,100]. Figure 6 shows the trajectory starting from this PRA 
extremum, whose initial coordinates are (1.690,5.380). For reference, two 
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Figure 7: Contours of the stream function (left), vorticity lj (middle) and 
potential vorticity q (right) of the modon solution (22). 


other trajectories are also shown with initial positions at 0.05 and 0.1 dis¬ 
tance from the vortex center. Due to the complexity of the flow, the trajec¬ 
tory patterns are not illuminating. However, their x- and ^-coordinates as a 
function of time reveal the oscillatory motion of nearby trajectories around 
the vortex center, while the vortex center itself has minimal oscillations (cf. 
middle and right panels of Fig. 6). The oscillations of the vortex center are 
due to the motion of the vortex as a whole. The nearby trajectories, how¬ 
ever, exhibit higher frequency oscillations which are due to their swirling 
motion around the vortex center. 

6.3 Stratified geophysical fiuid fiow 

We consider a simplified model for stratified geophysical fluid flow, the 
barotropic equation. This equation, in the vorticity-stream form, reads [41] 

dfCJ + + dx^ip = 0, cj = ( 21 ) 

where w{x^ t) and '0(x, 7/, t) are non-dimensional vorticity and stream func¬ 
tion, respectively. In deriving this equation, the viscous dissipation is ne¬ 
glected and the Coriolis frequency is assumed to be linear in the meridional 
coordinate y (i.e., the /3-plane approximation is used [41]). The Jacobian 
operator reads — dx^pdyUJ — dypjdxoo. The fluid velocity field u is 

given in terms of the stream function by u = {dypj^ —dx^ip)- 

Vorticity is not preserved along fluid trajectories when the flow satisfies 
(21). Instead, one can show that the potential vorticity q = uj-\-y is conserved 
along these trajectories (see, e.g., [41]). 
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Figure 8: The PRA for the modon solution (22). The initial time is to = 0 
and the final times are t = 100 (left), t = 250 (middle) and t = 500 (right). 
All figures are computed from a grid of roughly 35 thousand uniformly dis¬ 
tributed initial conditions in the unit disk. The points marked by crosses 
are the Lagrangian vortex centers obtained from Definition 3 over the cor¬ 
responding time intervals. 


We consider a steady exact solution of the barotropic equation (21) called 
a modon: a uniformly propagating vortex dipole. For this modon solution, 
the stream function and vorticity are given respectively by 


Hr, H 
Hr, v) 


(JH) \ . 
Ji{r) . 


0 < r < 1 

0 < r < 1 (22) 


in polar coordinates {r^Lp) where r = x‘^ + tan(/p = y/x and Ji is the 

Bessel function of the first kind [42]. This solution is written in a frame 
CO- moving with the modon at a constant speed c = 1. 

The stream function defines a flow on the invariant domain r < 1. 
While this solution can, in principle, be extended to the entire plane [42, 43], 
here we only consider the motion inside the unit disk. 

Figure 7 shows the stream function, the vorticity and the potential vortic¬ 
ity for the modon solution (22). Since the flow is integrable, its stream func¬ 
tion completely describes the flow structure, showing two counter-rotating 
vortices. 

The vorticity oj is negative in the upper half-disk y > 0 and positive 
in the lower half-disk y < 0. Its contours, however, do not reveal the two 
vortices present in the flow. This is because unlike the two-dimensional 
Euler flows, vorticity is not conserved along the trajectories of the solutions 
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of the barotropic equation (21). 

The potential vorticity g, as a conserved quantity, reveals the eddies. Its 
level curves (Fig. 7, right panel) resemble those of the streamlines. In fact, 
the particular solution (22) of the barotropic equation satisfies q — —'0. 

Figure 8 shows the PRA for integration times t = 100, 250 and 500. 
The integration time t — 100 is chosen such that almost all periodic orbits 
of X = u(x) complete at least one period. Even with this relatively short 
integration time, PRA contours already reveal the vortices. Obtained from 
a finite-time assessment of the flow, the PRA contours deviate from the 
trajectories. As the integration time increases, however, the PRA contours 
converge to the streamlines and Lagrangian vortex centers obtained from 
the PRA converge to the elliptic fixed points of the flow. 

6.4 ABC flow 

As our last example, we consider the Arnold-Beltrami-Childress (ABC) flow 
X = u(x) where 

( Asin( 2 ;) + Ccos(^)\ 

B sm{x)-\-Acos{z) j , (23) 

C sm{y) + B cos(x)/ 

with X = (x, z) and A, S, C G M are constant parameters [2]. The velocity 
field u is an exact steady solution of Euler’s equation for inviscid Newtonian 
fluids with periodic boundary conditions. The ABC velocity field is a Bel¬ 
trami vector field satisfying 6i;(x) = u(x) with cj = V x u being the vorticity 
field. 

In the following, we set A = 1, S = y^ 2/3 and C = ^1/3. The La¬ 
grangian computations are carried out on a uniform grid of 200 x 200 x 200 
initial conditions distributed over the domain G [0, 27r] x [0, 27r] x [0,27r]. 

Figure 9 (left panel) shows the helicity density (u,6t;) (=|6t; |2). While 
such Eulerian features may suggest coherent vortical motion throughout the 
domain, the ABC flow is known to have chaotic fluid trajectories in addition 
to coherent swirling trajectories lying on invariant tori [44]. 

These invariant tori form vortical regions that we seek to capture from 
finite-time flow samples as elliptic regions. Using a local variational principle 
extremizing Lagrangian shear, elliptic LCSs approximating the tori from 
finite-time flow samples have been constructed by Blazevski and Haller [14]. 
Here we illustrate that polar LCSs obtained from the PRA also give a close 
approximation at a reduced computational cost. 

Indeed, the PRA admits tubular level surfaces that closely approximate 
the invariant tori (Fig. 9, right panel). Codimension-two level sets of the 
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0 0 a: y 0 0 X 


Figure 9: Left: The helicity (u,6t;) = \uj\‘^ for the ABC flow. Right: The 
two-quadrant PRA 9l^ with the integration time t —to = 50, computed from 
formula (16). 


PRA are periodic material curves at the cores of the elliptic regions. These 
material lines serve as Lagrangian vortex centers by Definition 3. As in 
earlier examples, outside the elliptic islands formed by these closed level 
surfaces, PRA levels exhibit small-scale variations due to sensitive depen¬ 
dence of the rotation angle on the initial conditions. 

To examine how accurately the PRA field 9l^ captures the tori and the 
chaotic region boundaries, we release two trajectories from the initial con¬ 
ditions xo = (3.085,0,3.820) (red square in the bottom panel of Fig. 10) 
and Xo = (3.505, 0, 3.568) (blue square in the bottom panel of Fig. 10). The 
initial conditions are chosen such that they are nearby, yet one belongs to 
the chaotic region (red square) and the other (blue square) belongs to a 
smooth level-surface of the PRA signaling an invariant torus. 

These initial conditions are then advected under the ABC flow from time 
t = 0 to t = 500. The resulting trajectories are shown in the top panel of 
Fig. 10. As expected, the blue trajectory remains on a torus while the red 
trajectory exhibits chaotic behavior. Note that all curves correspond to a 
single trajectory and only appear as line segments because they are plotted 
modulo 27V. The intersections of the coherent trajectory with the Poincare 
section y = 0 shows that the PRA captures the invariant torus accurately 
(Fig. 10, bottom panel). 

We stress that both initial conditions studied here belong to topolog¬ 
ically equivalent regions of the local helicity (u,6t;) = The vorticity 

magnitude, therefore, fails to distinguish vortical regions from chaotic re- 
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Figure 10: Top: Two trajectories of the ABC flow. The blue trajectory 
starts in an elliptic island and traces the surface of an invariant torus. The 
red trajectory is chaotic. Note that the trajectories are plotted modulo 27r. 
The material line marking the Lagrangian vortex core (cf. Definition 3) is 
plotted in magenta color. Bottom: The initial condition (squares) of each 
trajectory is superimposed on the ^ = 0 slice of the helicity (u,6t;) = \uj\‘^ 
(left) and the PRA Ol^ obtained from formula (16) (right). The right panel 
also shows the intersections of the blue trajectory with the plane y = 0 (blue 
dots), as well as the vortex core (magenta cross). 


gions. This is because vorticity magnitude is not a material invariant of 
the Euler’s equation in three dimensions and therefore does not generally 
capture material behavior. 
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7 Conclusions 


Most approaches to coherent structures seek their signature in material sep¬ 
aration or stretching. By contrast, we have developed here an approach to 
locate coherent structures based on their signature in material rotation. To 
quantify finite material rotation in a mathematically precise fashion, we have 
used the polar rotation tensor from the unique rotation-stretch factorization 
of the deformation gradient. 

For two- and three-dimensional dynamical systems, we have derived ex¬ 
plicit formulas for the polar rotation angle (PRA) generated by the rotation 
tensor around its axis of rotation. While polar rotation has broadly been 
studied and used in continuum mechanics, the simple formulas we have de¬ 
rived here for the PRA in terms of the flow gradient, its singular values 
and singular vectors have not been available. These formulas enable the 
efficient computation of the PRA from basic quantities provided by existing 
numerical algorithms for Lagrangian coherent structure detection. 

Building on the PRA, we have also introduced the notion of polar La¬ 
grangian coherent structures (polar LCSs). These are tubular material 
surfaces along which trajectories admit the same PRA value over a finite 
time interval of interest. We have proposed regions filled by polar LCSs as 
rotation-based generalizations of the classic elliptic islands filled by KAM 
tori in Hamiltonian systems. 

As we demonstrated on a direct numerical simulation of two-dimensional 
turbulence, the PRA identifies Lagrangian vortex boundaries with high ac¬ 
curacy. While geodesic LCS theory of Haller and Beron-Vera [13] offers 
an exact detection of such vortex boundaries as solutions of differential 
equations, the present diagnostic detection of these boundaries as outer¬ 
most closed PRA level curves is substantially less computational, and hence 
preferable for an approximate identification of these boundaries. 

Outside the Lagrangian vortex boundaries, the PRA is dominated by 
small-scale noise due to its sensitive dependence on initial conditions. In 
these regions, therefore, the PRA displays no clear signature for hyperbolic 
LCSs governing chaotic tracer mixing. These latter types of LCSs, by con¬ 
trast, are efficiently revealed by another objective diagnostic, the finite-time 
Lyapunov exponent (FTLE) [15]. The PRA and FTLE have a well-defined 
duality: the former is a scalar field characterizing the rotational factor 
while the latter characterizes the stretch factor in the polar decomposi¬ 
tion of the deformation gradient. 

We have found that local extrema of the PRA mark initial positions of 
trajectories that serve as well-defined centers for elliptic islands. Oscillations 
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in these center trajectories are minimal and arise solely due to the material 
translation of the underlying island. Nearby trajectories inside the elliptic 
island, on the other hand, oscillate rapidly due to their swirling motion 
around the center trajectory (cf. Fig. 6). The ability of the PRA to identify a 
unique vortex center should be helpful in Lagrangian versions of the Eulerian 
eddy censuses carried out by Dong et al. [45] and Chelton et al. [46]. 

The elliptic island boundaries marked by PRA do not necessarily remain 
unfilamented under advection. If the goal is to find perfectly coherent La¬ 
grangian vortices (see, e.g., [47]), then the geodesic theory of Haller and 
Beron-Vera [13] should be applied. This theory identifies material vortex 
boundaries as closed null geodesics of the generalized Green-Lagrange strain 
tensor. The related computations require the a priori identification of phase 
space regions where such closed geodesics may exist [48]. Vortex regions 
identified from the PRA provide a quickly computable starting point for the 
detection of closed Green-Lagrange null geodesics. Incorporating the vortex 
centers obtained from the PRA in the geodesic LCS analysis is, therefore, 
expected to lead to a notable computational speed-up. 

Finally, the polar LCSs obtained as level curves of the PRA are frame- 
invariant for planar flows (see Proposition 3). Such objectivity is desirable 
for coherent structure identification methods in order to exclude false pos¬ 
itives and negatives specific to the coordinate system used in the analysis 
[15]. In three dimensions, however, the PRA does depend on the reference 
frame. The objective detection of higher-dimensional elliptic islands from 
their rotational coherence, therefore, requires further work. 


Appendix A Proof of Proposition 1 


Part (1): The trace of a tensor is independent of the choice of basis. If we 
represent the rotation tensor in the orthonormal basis {^k}i<k<3^ then 
its entires satisfy [R^q]^^- = Therefore, using formula (8), we 

can write 


‘r = E ((>’ !*■!.«<> = E («<. [UU ( 


i=l 

3 


i=l 






( 24 ) 


which, together with (12a), proves formula (13a). 
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To prove formula (13b), we first note the coordinate form of equation 

(12b): 


i (k: 




Ho] 


Applying the same argument used in (24) in the strain eigenbasis, we obtain 

*0 


sin 0*0 = 




Ho]/ 


* 7^ j- 


(25) 


Next we write the eigenvector in strain basis as fo obtain 


E = RJ. E [uy ■' E = E 

k k k k ^ 

which implies 


or, equivalently, = e, with and e defined in the statement of 

Proposition 1 . Since [^^tolk ~ formula (25) proves (13b). 

Part (2): Two-dimensional flows are parallel to a distinguished plane 
and exhibit no stretching or shrinking along the normal of this plane. In 
this case, we have 

Ai < A2 = 1 < A3, 

with the strain eigenvector ^2 pointing in the normal of the plane in question. 
Formula (13) then gives 


, <6,Vq„^3> 

cose^tg - - 1-7^-+ 


vV 


vV 


(26) 


Restricting our consideration to the two-dimensional plane of the flow, we 
reindex the quantities in formula (26) as A 3 ^ A 2 and ^3 ^ ^ 2 , given that 
the original A 3 strain eigenvalue of the flow is the second largest principal 
strain in the plane of the flow. After this re-indexing, equation (26) gives 


cos el 


= E 


1=1 




(27) 


The summands in this last expression are just the diagonal elements of the 
two-dimensional rotation tensor represented in the {^ 15 ^ 2 } basis (cf. 
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our discussion leading to equation (24)). Since the diagonal elements of any 
two-dimensional rotation matrix are equal, formula (14a) follows from (27). 
In two dimensions, the rotation tensor is of the form 


R*(xo)=(^ cos0^xo) sin0jg(xo)^ 


*0'' V - sin (xo) cos (xo) y ■ 

Thus, using the argument in (24), we obtain that 

sm< = (6,RS„6> = = 

whichi is the PRA formula (14b). 


(28) 


(6,VFUi 

a/m 


Appendix B Total Lagrangian rotation in planar 
flows 

The polar rotation 9\^ defined in Definition 1 is the net rotation of the 
{^ 1 ,^ 2 } eigenbasis over the time interval [to,t]. This quantity, however, mea¬ 
sures the rotation modulo 2tt and does not differentiate between rotation by 
00 and 00 + 2A:7r. Here, we also derive an expression for the total Lagrangian 
rotation of the eigenbasis that distinguishes between rotations differing by 
an integer multiple of 27r. 

Consider the equations of variations for a given infinitesimal displace¬ 
ment 

(29) 

Write ^(t) = reef) where = (cos 0, sin 0)^ and (r, 0) are functions of time. 
Substituting this in the equations of variations (29) we get 

re^ + r4)e-^ = rVu(x(f), t)©^, (30) 

with = (—sin</), COS0)"''. Since and are perpendicular, we have 

^ ^ {e^,Vu{x{t),t)e^), (31a) 

<i> ^ {ejyu{x{t),t)e^). (31b) 

Therefore, solving Eq. (31b), the total rotation of an arbitrary displacement 
vector ^0 = (cos 00 , sin 00 )^ is given by 

Otot ■■= 4>{t) -4>0= f Vu(x(r),r)e,^(^))dr. (32) 
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If the initial vector is chosen to be (or ^2)5 ^tot measures the total 
rotation of the eigenbasis {^1,^2}* We refer to 0tot as the total Lagrangian 
rotation. 

In practice, for evaluating the total Lagrangian rotation (32), one needs 
to first compute the deformation gradient from which the strain direc¬ 

tions {^15^2} are computed. The orientation of (or alternatively ^2) de¬ 
termines the appropriate initial condition = (cos 00, sin00with which 
Eq. (31b) should be solved. Note that Eq. (31b) must be solved simulta¬ 
neous with the dynamical system x = u(x,t) since Vu is evaluated along 
trajectories x(t;to,xo). 

Therefore, evaluating the total Lagrangian rotation is more expensive 
than computing the PRA. The connected components of the level sets of 
Oiot and are identical by an argument similar to the one used in the 
proof of Proposition 2. Thus the polar LCSs revealed by these two scalars 
are also identical. 


Appendix C Proof of Proposition 3 


Differentiating both sides of the formula (17) with respect to the initial 
condition xq gives 

= Qit)VFlQ^{to), (33) 

where denotes the deformation gradient in the y = K (yo) coordinate 

system. From (33), we obtain 

VF;„ = QT(t)VF*„Q(fo) 

= QT(t)R*^Q(to)QT(to)U*„Q(to) 

= 

where the rotation tensor ^^d the positive definite, 

symmetric tensor {to)\Jl^Q{to) represent the unique polar decom¬ 

position of VF^q . Then 


tr^Uyo) = tr 


Q"^WR'Uxo)Q(to) 


^ f cos [el (xo) + q{to) - q{t)] - sin [Ol (xq) + q{to) - q{t)] \ 

\sm [(9^xo) + q{to) - q{t)] cos [6*^x0) + q{to) - q{t)] ) 

= 2 cos [6*40 (xo) + qito) - q{t)] , 
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where q{t) represents the angle of rotation associated with Q(t). Therefore, 
if the polar rotation angle generated by transformed rotation tensor is 
Otoiyo), then 


cos (0to(yo)) = itrR^yo) 

= cos {9l^{xo) + qito) - qit)) . (34) 

Consequently, if two points xq and xq lie on the same connected level 
set of 0 ^q(xo), then the corresponding points also lie on a connected level set 
of 0fo(yo), even thongh we generally have ^ dj^iyo)- 

We note that the level sets of PRA in three dimensions are generally not 
objective. An essential part of the above argument, leading to equation (34), 
is that the rotation matrices Q(t), R^q(xo) and Q(to) share the same axis 
of rotation (i.e., the normal to the plane of motion). In three dimensions, 
such a uniform axis of rotation does not generally exist, and hence a relation 
similar to (34) does not hold. 
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